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ABSTRACT 

We have developed a new halo finding method, Physically Self-Bound (PSB) group finding algorithm, 
which can efficiently identify halos located even at crowded regions. This method combines two 
physical criteria such as the tidal radius of a halo and the total energy of each particle to find member 
particles. No subtle dependence of halo mass functions on various parameters of the PSB method 
has been found. Two hierarchical meshes are used to increase the speed and the power of halo 
identification in the parallel computing environments. First, a coarse mesh with cell size equal to 
the mean particle separation i mea n is used to obtain the density field over the whole simulation box. 
Mesh cells having density contrast higher than a local cutoff threshold <5loc are extracted and linked 
together for those adjacent to each other. This produces local-cell groups. We analyze the group of 
particles located at each local-cell group region separately. This treatment makes the halo finding 
method easily implemented on the parallel computing environments since each computational rank 
takes the halo identification job in each particle group independently. Second, a finer mesh is used 
to obtain density field within each local-cell group and to identify halos. We set the cell size of the 
refined mesh to be twice the gravitational force softening length e. The density peaks in the fine 
mesh are the halo candidates. Based on the fine mesh, we split particles into many groups located 
in the density shells with different density levels. If a density shell contains only one density peak, 
its particles are assigned to the density peak. But in the case of a density shell surrounding at least 
two density peaks, we use both the tidal radii of halo candidates enclosed by the shell and the total 
energy criterion to find physically bound particles with respect to each halo. 

We have tested the PSB using a binary halo model against other popular halo finding methods, such as 
the Friend-of-Friend (FoF), Spherical Overdensity (SO), DENMAX, and HOP. Similar to DENMAX 
and HOP, the PSB method can efficiently identify small halos embedded in a large halo, while the 
FoF and the SO do not resolve such small halos. We apply our new halo finding method to a 1- 
Giga particle simulation of the ACDM model and compare the resulting mass function with those of 
previous studies. The abundance of physically self-bound halos is larger at the low mass scale and 
smaller at the high mass scale than proposed by the Jenkins et al. (2001) who used the FoF and SO 
methods. 

Subject headings: Cosmology: N-body simulation: halo : halo-finding methods: Numerical 



1. INTRODUCTION 

Cosmological N-body simulations (Park 1990; Gelb & 
Bertschinger 1994; Park 1997; Evrard et al. 2002; Du- 
binksi et al. 2003; Bode & Ostriker 2003, among others) 
have been used to test cosmological models in various 
fields of interests. To compare simulation results with 
observed galaxies or other visual cosmic objects, one has 
to extract virialized dark matter halos from the distribu- 
tion of simulation particles. 

An easy way to identify virialized halos from particle 
distribution is to link particles with distances less than 
^FoF (Friend-of-Friend, Audit et al. 1998; Davis et al. 
1985). The value of the linking length ipoF is usually 
set to 0.2 x /mean corresponding to the virialization over- 
density p/pi, — 178 (Porciani et al. 2002), where pi, is 
the mean background density. The chain of linked parti- 
cles forms a group of particles and the particle group is 
considered as a virialized halo. But the FoF cannot iden- 
tify small halos embedded in larger high-density regions. 
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This effect is similar to the "overmerging" problem oc- 
curred in the poor-resolution simulations (Moore et al. 
1996; Klypin et al. 1999). The main drawback of the FoF 
is the "overlinking" (Gelb & Bertschinger 1994) since, in 
a binary or multiple halo system, member halos are often 
linked by bridging particles. Then a dumbbell-like halo 
is identified and, consequently, the halo quantities, such 
as the center of mass, the shape, the rotation velocity, 
etc., are blended. The Spherical Overdensity (SO, Lacey 
& Cole 1994; Warren et al. 1992) uses the mean density 
for virialized halos. The SO searches for density peaks 
and puts spheres around them by increasing the radius 
of the sphere until the internal mean density satisfies the 
virialization criterion (p/pb = 178). Particles inside the 
sphere are grouped as members of the spherical halo. 
However, small halos embedded in virialized regions can 
not be properly resolved. 

More improved halo finding methods have been pro- 
posed in the last decade. The particle sliding scheme 
in the density field is adopted by the DENMAX (Gelb 
& Bertschinger 1994) or its variant SKID (Weinberg et 
al. 1997; Jang-Condell & Hernquist 2001), and the HOP 
(Eisenstein & Hut 1998). The DENMAX uses rectan- 
gular density grid cells to slide particles to the nearby 
densest grid cell. Then it scoops up particles that are 
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stacked at local density maxima and checks the total en- 
ergy of each particle to discard unbound particles. The 
SKID method, an improved version of the DENMAX, 
uses a variable smoothing length and the density gradi- 
ent in a coordinate-free density field. The HOP calcu- 
lates the density field in the way similar to the SKID. 
However, it uses a different type of particle sliding. The 
HOP method searches for the maximum density among 
a particle's nearest neighbors. Particles are slid into the 
nearby densest particle. The HOP groups particles in 
local density maxima as virialized halos similar to the 
DENMAX. 

But above halo finding methods arc mainly based only 
on the density field in determining the member particles 
of halos. Namely, used are the density-related quanti- 
ties such as the linking length (FoF), the overdensity 
(SO), and the sliding on a local density field (DENMAX, 
SKID, and HOP). As a result, the distribution of parti- 
cles in a halo depends on the density geometry. For ex- 
ample, the FoF is likely to produce dumbbell-like halos 
and the SO finds halos with spherical boundaries. Both 
methods cannot resolve any substructure inside the viri- 
alized region. Using the particle sliding scheme seems 
more proper than the overdensity criteria (p/ pb and Zfof) 
since it is reasonable to assume that boundaries of viri- 
alized halos may follow the 3-dimensional density valley 
when they are next to other halos. But, in some cases, 
the density gradient field does not properly describe halo 
boundaries when there is no density valley. For example, 
consider a small halo is located in a large halo. Around 
the small halo, there may exist a density valley between 
the two density peaks of halos. However, there is no 
density valley in the outskirts of the large halo region 
beyond the small halo. So the density valley does not 
form a closed surface. Consequently, particles that are 
members of the large halo may slide to the small halo 
along the path without encountering any density valley. 

In this paper, we present a new halo finding method 
based on the density map. However, our method also 
takes into account the tidal radius that forms a closed 
surface around the halo and the total energy to select 
physically bound particles only. We describe how to di- 
vide simulation particles into local particle groups in §21 
For individual particle groups, we present the way to 
identify virialized and stable halos in the tidal field in 
fJ3] The values of the group finding parameters are given 
in 21 and we test the sensitivity of the PSB to param- 
eters in <J3 In fjni we compare our halo finding method 
with other popular methods. In §T| we apply our method 
to the I-Giga particle simulations. Conclusions follow in 

m 

2. LOCAL PARTICLE GROUPS ON A COARSE MESH 

The PSB is based on the two hierarchical meshes for 
particle-density field. A low-resolution density field is 
used to find the "local particle groups" , and for each lo- 
cal particle group we use a high-resolution density field 
to find density peaks and to divide particles into multiple 
"particle sets". The division of whole simulation parti- 
cles into many local particle groups allows us to reduce 
computational costs in halo findings. 

The PSB uses all simulation particles to build the den- 
sity field on a coarse mesh with cell size equal to the 
mean inter-particle separation l mca , n . To assign densities 



to the mesh cells we use the W& (Monaghan & Lattanzio 
1985) 

i f (1- §z + jx 3 ) if < x < 1, 
Wi(r,h) = - I \{2~xf ifl<ar<2, (1) 

[ otherwise, 

where r is the distance of a particle to a grid point, h is 
the smoothing length, and x = r/h. We set h equal to 
Wan in the coarse mesh. We attach particles to mesh 
cells as the linked list (Hockney & Eastwood 1981). We 
extract overdense cells enclosed by an isodensity surface 
defined by the density contrast threshold <5loc- These 
cells are grouped and labeled. To these grouped cells, 
we attach adjacent underdense mesh cells to avoid sharp 
boundary truncations of halos. The particles in a local- 
cell group are found and called the local particle group. 

3. PARTICLE SET HIERARCHY AND PARTICLE 
ASSIGNMENT 

3.1. Local Density Field on a Fine Mesh 

For each local particle group, we construct a high 
resolution density field on a fine mesh with cell size 
and smoothing length equal to twice the force soften- 
ing length e to resolve small but physically meaningful 
structures. We search for local density peaks above the 
density contrast threshold 8 p . Those peaks are called 
halo cores. Around each density peak, we find the min- 
imum density level S c with which the isodensity surface 
encloses only that peak (see Fig. 1). Particles in cells 
surrounded by this isodensity surface are set to be core 
members of the density peak. We constrain that the 
number of core particles (Af core ) should be at least 10 
to avoid the detection of spurious peaks induced by the 
Poisson fluctuation. We assume those core particles as 
members of the halo. Other remaining particles are di- 
vided into subsets located within the density shells hav- 
ing Ni evel density levels between <5 LO c and S CjTnax , the 
highest S c . Every density shell encloses at least two halo 
cores. We call the group of particles in a density shell 
a particle set (PS). Now we have to assign each of these 
particles to one of the halo cores. 

We illustrate our halo finding method in Figure ^ 
Shown are two regions surrounded by the outermost den- 
sity contours with the threshold level (5loc • As described 
in the previous section we call each region a local-cell 
group. The boundaries of the inner-most density region, 
PI, P2, and P3 are found by searching for <5 c 's. Parti- 
cles in those regions are the core members of each density 
peak. Now we call each of PI, P2, and P3 a halo can- 
didate. To simplify the illustration in Figure 1., we use 
three density threshold (Ni eve i — 3) and get three parti- 
cle sets spatially split as Al, A2, and A3. But actually, 
we use Ni eve i = 10- Particles in the region Al may be 
assigned to one of PI and P2. Particles in the regions 
of A2 and A3 can be, similarly, assigned to one of the 
three halos (PI, P2, and P3). 

3.2. Tidal Radius Criterion 

The PSB method uses multiple steps to assign par- 
ticles to halos. Particle sets are sorted in a decreasing 
order of density levels. Starting from particles of the top 
inner-most PS, we apply the total energy check to each 
particle and the tidal radius criterion to each halo en- 
closed by the PS. But there needs an iterative scheme in 
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Fig. 1. — Schematic illustration of a density map in our halo 
finding algorithm. In panel (a), The regions PI, P2, P3 and A4 
contain only one density peak with <5 > 8 P , and are enclosed by 
isodensity contour surfaces. The regions A3 and A4 are identified 
by the density threshold <5loc- In panel (b), we show the density 
map along the line x in panel (a). A horizontal dotted line shows 
the core density contrast <5 C of each peak. We use N[ eve [ = 3 and 
show the density level in a horizontal dashed line. 



this hierarchical membership determination. Since a halo 
mass is usually measured by its member particles, a halo 
potential may be underestimated unless all PS's down to 
the bottom one are taken into account. To supplement 
the halo mass, we use the PS's having the density levels 
both higher and lower than that of the current particle 
set PS C . Let the PS; be the particle sets geometrically 
enclosing the halo and having density levels lower than 
that of the PS C . The potential Q(k) of the fc'th particle 
with respect to its host halo is calculated by 

*(* e PSc) = - £ Gm,_ £ ^, (2) 



where r k i — \ fe~ f%\, rri£ is the tth particle mass, MEM 
are the member particles already assigned to the halo 
in the PS-f's which have density thresholds higher than 
that of the PS C . The second term of the right hand side 
of equation J2J) significantly contributes to the particle 
potential $ if the number of particles in PSj's is large. 
We do not include particles of the PS C for the halo mass 
in equation J3J) since an accidentally close particle to the 
fc'th particle can make the potential overestimated. 

In some cases, particles happen to be bound to many 
halos simultaneously. To resolve this situation, we pro- 
pose a new halo boundary indicator, the tidal radius, 
which is based on the local gravitational field. A tidal 
radius of a halo is calculated from a tidal mass M t which 
is the number of only the halo member particles 



M t 
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keMEM 



(3) 



Let R t (i,j) be the tidal radius Rt(i) of the i'th halo with 
respect to the j'th halo more massive than the i'th. A 
simple analytic form for the tidal radius is the Jacobi 
radius Rj 



Rj(i,j) = D 



M t {i) 



3M t (j,r < D) 



1/3 



(4) 



where D is the distance between the centers of mass of 
the halos, and M t (j, r < D) is the tidal mass of the j'th 
halo mass contained in a sphere with a radius D. The 
effect of underestimation of the halo mass (Eq. EJ) on the 
tidal radius (Eq. ^} is thought to be insignificant due 
to three reasons. First, Rj is proportional to its cubic 
root. Secondly, we measure M t at each PS analysis and 
iteratively update the tidal radius of each halo. Then the 
tidal radius of a halo converges to a true value when the 
bottom PS is considered. Thirdly, the j'th halo member 
particles beyond D do not contribute both to M t (j,r < 
D) and M t (i). The particles of a PS; that surrounds 
both the j'th and j'th halos are usually too far from 
the two halos to contribute to the tidal masses. And it 
is less likely for the particles of the PS; to be assigned 
to the i'th halo and, M t (i) is not expected to increase 
significantly. Therefore, even if we initially assign an 
inaccurate value to the tidal radius, the tidal radius of 
each halo converges to a true value in this hierarchical 
particle assignment process. 

However, the Jacobi radius Rj (Eq. is obtained 
assuming that M t (i)/M t (j,r < D) « 1. To take into 
account more general situations, we calculate the tidal 
radius Rt of a smaller halo as a function of the halo mass 
ratio m = M t (i)/M t (j,r < D) using 



(1 - r) - mr" 2 - 1 + (1 + m) r = 0, 



(5) 



where v = R t /D (Eq. 7-82 of Binney & Tremaine 1994). 
We assume that two halos orbit circularly around their 
center of mass. In Figure [21 we show the ratio of the 
tidal radius Rt to the distance D as a function of the halo 
mass ratio m under various conditions such as the circu- 
lar orbital motion (Eq. [5J solid curve), the Jacobi limit 
m — > (dotted curve, King 1962; Binney & Tremaine 
1994), and no orbital motion (short-dashed curve). As 
can be seen in Figure our adopted tidal radius model 
(solid curve) can be approximated by the Jacobi radius 
Rj in the limit of m — » 0. However, as m — + 1 in a binary 
system with equal halo masses, the Jacobi radius Rj ap- 
proaches an incorrect value of 0.69Z3 while our formula 
gives Rt = D/2. We fit our model to a fitting function 
(Keenan 1981) 



Rt = D 



M t {i) 



a(M t (i) + M t (j,r < £>)) 



(0 



where a = 4.813 and j3 = 0.318 (long-dashed curve). 

To a small halo many tidal radii can be assigned with 
respect to other more massive halos enclosed by a given 
PS. We adopt the minimum value of i?t(i)'s as 

R t {i) = MLN{R t (i,j = l),.-.,R t (i,j = N eh )\, (7) 

where N e h is the number of enclosed halos with mass 
satisfying M t {j,r < Dij) > M t (i), where Dy is the dis- 
tance between the i'th and j'th halos. If a particles is 
not bound to any halo or beyond the tidal radius of each 
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m=M t (i)/M t (j,r<D) 

Fig. 2. — r-m relation. We show the ratio t of the tidal radius to 
the distance between halos as a function of the mass ratio between 
Mt(i) and Mt(j, r < D) for the circular orbit (solid curve), for the 
Jacobi limit (dotted curve, Eq. |1J, and for the system with no 
orbital motion (short-dashed curve). The long-dashed curve shows 
the fitting result of the tidal radius for a circular orbit system 
adopted in the PSB method. 



halo, we stack it in the free-particle list. Particles in the 
free-particle list are temporarily added to the particle list 
of the next PS. 

After completing the above steps for all PS's we use the 
Friend-of-Friend (FoF) method with ZpoF = 0.2 x Z mcan 
to find particles belonging to virialized halos. Particles 
which are not members of any halo are stacked in the 
free-particle list again. 

Before we get final halo data, we check whether or not 
halos are self-bound and virialized. This process is simi- 
lar to those of other halo finding methods except for the 
tidal radius criterion. We find the minimum tidal radius 
of each halo using equation (JJJ. At this time, N e h is the 
number of all halos more massive than a given halo in the 
local-cell group region. We look into the member-particle 
and free-particle lists. Unbound particles or those be- 
yond the tidal radius of a halo are discarded and stacked 
to the free-particle list. We iterate this process four times 
to obtain self-bound and stable halos. Finally, the FoF 
method is applied to satisfy the virialization condition of 
each halo. 



The PSB method is summarized in the flowchart in Fig- 
urc|21 In panel (A), we show the pre- halo finding process 
described in SJ21 T° find local particle groups, we calcu- 
late the density on a low-resolution mesh using all simu- 
lation particles. Using densities on the mesh and a pre- 
defined parameter <5loCj we find overdense cells and link 
them making local-cell groups. Particles on the linked 
cells are grouped and extracted forming local particle 
groups. Using a local particle group, the PSB method 
identifies virialized halos. Panel (B) shows a flowchart 
of the halo finding process on a fine mesh. After find- 
ing particle sets, we iteratively update the tidal radius 
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Fig. 3. — Flowchart of the PSB method. We show the way 
to extract local particle groups based on a coarse mesh described 
in panel (A), and to identify virialized halos from a local particle 
group in panel (B). There are four critical parameters used in the 
PSB method, such as <5loCj ^p> N CO re, and Ni eve l- 



of a halo and calculate the total energy of member par- 
ticles. We, then, apply the FoF method to find member 
particles which satisfy the virialization condition. 

3.3. Parallelization 

We parallelize our halo finding method using the 
Message-Passing-Interface (MPI) library. As described 
previously, the PSB method is based on two meshes with 
coarse and fine cells. It therefore adopts two different 
modes of domain decomposition. Since a coarse mesh is 
built over the entire simulation box, it is desirable for 
each computational rank to have the same domain size. 
We use the domain decomposition into the z-directional 
domain slabs as described in the GOTPM code (Dubinski 
et al. 2004). We search for overdense cells and connect 
them together. Overdense cells at the upper and lower 
boundaries of local domains can be linked to neighboring 
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domains. Because we use one-dimensional cyclic order- 
ing of local domains, we pad images of mesh slices im- 
ported from neighboring domains at the upper and lower 
boundaries of the local domain. Those imported mesh 
slices contain local informations of their domains. We 
also take account of the periodic boundary conditions of 
the domain slab in the horizontal directions. With these 
linked groups of cells, we pick up particles in the over- 
dense cell regions and save those particles to disk. 

On the fine mesh the domain decomposition is made 
based on the local particle group. We set a master rank 
to do the disk I/O of the local particle data and to dis- 
tribute them into one of the slave rank that is idle at 
that time. Each slave rank receives particle data and 
performs the process described in the JO] and H3.2I Af- 
ter finishing the halo finding procedure in a local particle 
group, the slave rank sends halo member particles back 
to the master rank to save the halo finding results. 

4. PARAMETER DETERMINATION 

The lower Shoe is, the smaller are the smallest halos 
identified. Two important factors constraining Shoe are 
the computer memory available for each rank and the 
maximum allowable number of density peaks in a local 
particle group. Since a fine mesh devours much memory 
budget, the mesh with about 750 3 cells is the maximum 
allowable size on 32bit machines. Second, the tidal radii 
should be measured for a larger number of density peaks 
when there are many density peaks in a local particle 
group. Counting particles to measure M t (j,r < D) in 
eq. © with respect to every smaller halo consumes a lot 
of CPU time. 

We have empirically found that it is adequate to set 
Shoe = 1-3(72(0), namely, ^loc ~ cr 2 (0) for the case 
of the SCDM and <5loc ~ 3(72(0) for the case of the 
concordance ACDM model where (72(0) is the standard 
deviation of density fluctuation measured in the coarse 
mesh at z = 0. We use a constant value of Shoe at 
all redshifts. A lower value of <5loc enable one to iden- 
tify smaller halos in the lower density regions, but the 
volume of each local-cell group increases significantly. In 
the extreme case of <5loc = — 1, there is one big local-cell 
group containing all simulation particles. Therefore, the 
advantage of dividing all simulation particles into mul- 
tiple local particle groups is no more sustained. On the 
other hand, higher value of Shoe is not desirable since the 
tidal radii can be underestimated due to the insufficient 
description of the local tidal field. Additionally, small ha- 
los residing in lower dense regions can be missed in the 
process of the halo findings. We choose larger <5loc f° r 
the ACDM model because the ACDM model has higher 
degree of connectivity of overdense regions (filamentary 
structures) than the SCDM model. Above values of Slog 
are determined for 2GB memory budget. If more mem- 
ory budget is available, we are able to reduce Shoe and 
identify smaller halos in lower dense regions. 

To avoid detecting spurious density peaks, we require 
density peaks to have at least two particles in the fine 
mesh cell. The density corresponding to this constraint 
is 

^2x(Un///) 3 , (8) 

where If is the cell size of the fine mesh. If l mean /lf = 5, 
then S ~ 250. A halo having the density peak cell with 
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two particle satisfies the virial condition S v 
matically. 

satisfy the virial condition 

^LOC &vir < S p . (9) 

In the PSB method, two density meshes with different 
cell sizes are adopted and different density thresholds, 
Shoe and S p , are applied to them. Occasionally, there 
may be peaks with S > 5 P at the scale of the fine cell 
but with 5 < Slog on the coarse cell. We estimate the 
fraction of such missing peaks in the following to show 
the problem is minor for our choice of <5loc and S p . 

For a Gaussian density field, the probability of finding 
the regions with S > 5 t is 

p[iai ' 1 = ^f * exp (£)' (10) 

where a 2 is the variance of the density fluctuation. Let 
us use the subscript 1 to denote the statistical quantities 
measured in the fine mesh and the subscript 2 in the 
coarse mesh. In our halo finding method the smoothing 
lengths h's are set to the mesh cell sizes, namely hi — 



li = 2e and /i 2 



h 



It is usually e = 0.1 x 



1, 



in high resolution cosmological simulations. The 



conditional probability that a region has S > S p in the 
fine mesh when it has S < <5loc in the coarse mesh is, 
P[S hl >S P \S h2 <S hOC ] 

= P [S hl > S p , 6 h2 < Shoe] IP [S h2 < W] • (11) 
The joint probability P [Sh ± = Si,Sf l2 = S 2 ] is 

p [Shi = 8i>8h 3 = <y 
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(13) 
(14) 



a =(7, - 



d 3 kP(k)W^(k)W hl (k), (15) 



(16) 



and P(k) is the power spectrum of the density field. The 
smoothing kernel Wk in fc-space is given by 
1 f°° 

W h (k) = — - / dr 3 W 4 (r- h) exp (tk • r), (17) 



(2tt) 3 Jo 

and we use a spherically symmetric smoothing window 
Wa (Monaghan & Lattanzio 1985). Then, we get 



P[Shi >S P \S h2 <S L oc] 
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(18) 



6 



Kim and Park 



where 

N = P[S h2 <W] 

=s/, *"(-59)- (19) 

We have applied our halo finding algorithm to two 
cosmological N-body simulation data. The simulations 
are based on a ACDM and the SCDM models with the 
same simulation box size Li,ox = 512/i _1 Mpc. In the 
ACDM and the SCDM models we follow the gravita- 
tional evolution of 1024 3 particles in 1024 3 mesh cells 
from z — 23 to z = with 460 time steps by using the 
GOTPM code (Dubinski et al. 2004). We set S LO c = 10 
(~ 3ct 2 (0)) for the ACDM model and S LOC = 5 (~ cr 2 (0)) 
for the SCDM model, and fix the peak density thresh- 
old dp = 312.5. We calculate the number of unidentified 
density cells (<5i > 5 P ) which reside in the underdense 
background (<j 2 < Slog) using equation (JTSJ) . At z = 
the number of such density peaks on the fine mesh in 
the total simulation box is only a few, which is very 
small compared to the number of halos (~ 2 x 10 6 in 
the ACDM simulation and - 4 x 10 6 in the SCDM sim- 
ulation) identified by the halo finding method. 

In some situations, small halos identified in overdense 
backgrounds may have the true peak densities lower than 
S p when the underlying backgrounds are removed. To 
obtain unbiased halo samples we calculate the halo den- 
sity by using their member particles at the end of halo 
identification, and discard the halos having the peak 
density S < 5 P from the halo catalog. Finally, for our 
ACDM and SCDM simulation data we obtain halos with 
masses larger than 6 x 1O 11 /i _1 M0 (~ 58 particles) and 
10 12 /i _1 M© (~ 30 particles) , respectively. 

5. SENSITIVITY TO PARAMETERS 

Using halo mass functions, we examine the dependence 
of halo finding results on the four group finding param- 
eters, <5 L oc, S' p , Ni eve i, and N core , where S' p = 5 p /125. 
As a testbed, a ACDM model is simulated with 128 3 
particles in a box of size 128/i~ 1 Mpc. The bias fac- 
tor is set to b = 1.11 and the initial epoch is Zi = 4. 
At z = 0, the RMS of linear density fluctuation on 
the coarse mesh is 2.69. We use the canonical param- 
eter set (S LO c,S'p,N leve i,N core ) = (11,2.5,10,10). Then 
S' p ~ 4.1 x tr a (0). 

In panel (a) of Figure 01 we show halo mass functions 
with varying S' p from 1.5 to 2.5. When lower values of 
Sp are adopted, much more smaller halos are found. In 
panel (b) the halo mass function shows negligible depen- 
dence on N core . As the peak parameters, S p and N core 
are correlated with each other. The variation of N core 
apparently takes no effect on the halo mass functions 
since S p is a more stringent parameter when S p is high 
enough to satisfy 5 P > S V i r . The dependence of the halo 
mass function on the threshold density 5loc is also weak 
as shown in panel (c). The parameter Ni eve i hardly af- 
fects the halo mass functions as demonstrated in panel 
(d). We use a high value of Ni eve i = 10 to better resolve 
clouded regions with many peaks. For example, there are 
about two thousands peaks in the largest local particle 
group in our 1024 3 particle simulations at z = 0. 

We conclude that the most important parameter is S p 
which determines the lower limit of mass of the halos 
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Fig. 4. — Halo mass functions with various different 
parameter sets. The canonical parameter set is obtained 
(5 hoc ,6' p ,Ni eVBl , Ncore) = (11,2.5,10,10), where 5' p = <5 P /125. In 
panel (a), we show the halo mass functions varying 8' p from 1.5 to 
2.5. Panel (b) is same as panel (a), but N core is varied. We vary 
<5loc m P ane l ( c ) an d Ni eve i in panel (d). The solid line shows the 
Sheth Sz Tormen mass function, and the two vertical arrows mark 
the mass of 50 and 100 particles. 



found. The role of <5loc m the PSB method is less sig- 
nificant when S p is high. But if one wants to identify less 
massive halos or to reduce S p , attention must be paid 
to the value of <5loc- He must measure the probability 
of missing fine cells which have 5 > S p in the under- 
dense background (S < #loc)- The other parameters are 
proved to be insensitive to the halo finding results in the 
PSB method. 

6. COMPARISON WITH OTHER HALO FINDING 
METHODS 

We have compared the results of the PSB method with 
those of other competing halo finding methods. We con- 
struct a binary system of halos having 20000 and 1000 
member particles, respectively. The large halo has the 
virial radius R — (SMh/^irpng) 1 ' 3 , and the small halo 
has its tidal radius under the tidal field of the large halo. 
The distribution of the member particles is set to follow 
the isothermal density profile. The small halo is slightly 
offset in x-direction (upper- left panel in Fig. [SJ. The 
directions of particle velocities are chosen randomly and 
the magnitudes are given in accordance with the viri- 
alization condition. Here, the small halo is bound to 
the potential well of the large halo, with no bulk mo- 
tion. The upper-left panel of Figure shows the particle 
distribution of each halo with the enlarged particle dis- 
tribution of the small halo in the inset. We also mark 
the small halo region with a box in the panel but do not 
overplot the member particles of the small halo for clar- 
ity. The middle-left panel shows the result of the FoF 
method which does not identify the small halo. The halo 
boundary found by the FoF method is not smooth due 
to the Poisson fluctuation of particles at the outskirts 
of the large halo. The SO method cannot identify the 
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Fig. 5. — Comparison among the results from various halo find- 
ing methods. We show the model distribution of particles of two 
halos in the upper-left panel. The large halo has 20000 particles 
following the isothermal density profile. The small halo whose re- 
gion is marked by a box has 1000 particles with the isothermal 
but tidally truncated density slope. In the inset, we show the en- 
larged view of the small halo for clarity. The small halo is offset to 
right with respect to the large halo. We show the identified parti- 
cle distribution obtained by halo finding methods, such as the FoF 
(middle-left panel), SO (middle-right panel), DENMAX (lower-left 
panel), HOP (lower-right panel), and PSB method (upper-right 
panel). The insets at the upper-right corners show the small halos 
identified by these methods. 



small halo like the FoF method. The large halo is found 
with a spherical boundary, which is a characteristic of 
the SO method. These two methods cannot detect the 
small halo since the halo identification is mainly based 
on the mean overdensity criterion. 

The lower-left panel shows the result of the DENMAX 
method, which resolves the small halo. However, due 
to the sliding of particles to the nearby densest cell the 
large halo seems to have no member particles at the right 
hand side of the small halo region. All those particles 
are assigned to the small halo, but they are not bound 
because of the shallow potential well of the small halo. 
Consequently, those particles are regarded unbound, and 
erased from the member list of the small halo. Having 
the particle sliding algorithm similar to that of the DEN- 
MAX, the HOP method can detect particles beyond the 




30 

xi.. -M t .< : ; 

Fig. 6. — Projected particle distribution in the ACDM simu- 
lation. The thickness of the slice is 10h~ 1 Mpc. Only every tenth 
particles in the volume are shown for clarity. Projected ellipsoidal 
shapes (marked by ellipses) of total 355 halos identified in this 
volume arc shown. 



small halo region. Those particles can slide according to 
local particle links, and possibly make a detouring path 
around the small halo region to the density maxima of 
the large halo. Since we do not use the REGROUP algo- 
rithm to regroup particles stacked in the local and small 
density maxima, there are many "particle holes" in the 
large halo identified by the HOP result. 

The upper-right panel shows the result of our PSB 
method. The boundary of the large halo is similar to that 
in the FoF method because the FoF method is used in our 
algorithm to find member particles in the virialized re- 
gion. The PSB method identifies member particles of the 
small halo nearly the same as the input. For examples, 
the numbers of identified member particles of the small 
halo are 0, 0, 758, 1101, and 1091 for the FoF, SO, DEN- 
MAX, HOP, and PSB methods , respectively. Therefore, 
the PSB method seems to recover the input particle dis- 
tributions more accurately than other methods except for 
the HOP. But the density allocation method of the HOP 
requires more CPU time to search for N ne i g hbor parti- 
cles than the PSB method which uses a fixed smoothing 
length. We conclude that the PSB method is an efficient 
and accurate halo finding method that can be adapted 
to the parallel computing environments. 



7. RESULT OF APPLICATION 



7.1. 



Halo Finding for 1-Giga particle simulations 

We apply our halo finding method to the 1-Giga par- 
ticle data simulated on the IBM SP3 at Korea Institute 
of Science and Technology Information (KISTI) and on 
a Linux cluster. The times needed to complete the halo 
findings at z = in the SCDM and the ACDM mod- 
els are about four hours when 32 POWER-4 CPUs are 
used. In Figure El we show the result of halo finding in 
the ACDM model at z = 0. Only every tenth particles 
in the volume are shown for clarity. The projected ellip- 
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Fig. 7. — A mock redshift survey in the ACDM model mimicking 
the LCRS. The simulated survey slice is 360° by 4.5° wide. Total 
number of mock galaxies is 27958, with mass 
The highest redshift z max is 0.02. 



soidal shapes of halos are plotted on top of the particle 
distribution. 

7.2. Mock Galaxy Redshift Survey 

Figure \7\ shows the distribution of halos in a mock 
galaxy redshift survey similar to the Las Campanas Red- 
shift Survey (hereafter LCRS, Lin et al. 1996). The 
selection function adopted here is the same as that of 
the LCRS. We do not take into account the conversion 
of halo mass functions to galaxy luminosity functions 
(Yang et al. 2003). The effects of peculiar velocities 
of halos are included. The halo number density in our 
ACDM simulation is p h = 0.0082/i 3 /Mpc 3 which is about 
3.5 times less than that (p g = 0.029 ± 0.002/i 3 /Mpc 3 ) of 
the LCRS galaxies in absolute magnitude range —23.0 < 
Ai — 5\ogh < —17.5. Using the Schechter function 



(f>(M) = (QAhil0)<t>* 



-^q0A(M* -M) 



(1- 



x exp 



^q0A(M* -M) 



(20) 



where M* = -20.29, a = -0.7, and (jy* = 0.019 (Lin 
et al. 1996), we find that the mass 6 x 10 11 /i _1 M Q cor- 
responds to the magnitude —19.8 + 51og/i in the LCRS 
galaxies from a direct comparison between the cumula- 
tive number density of the LCRS galaxies and that of 
the ACDM halos. 

7.3. Isolated and Distinguished Halos; Halo Mass 
Function 

In this subsection, we compare the halo mass func- 
tions found by the PSB and the FoF-like methods with 
the well-known analytic and fitting functions. We call 
the halo found by the PSB the distinguished halo. A 
distinguish halo is a virialized structure irrespective of 
its environment. For example, each halo in Figure [B] is 



regarded as a distinguished halo. The FoF-like halo cat- 
alog is constructed by merging halos when halo regions 
are overlapped. To model each halo as an ellipsoid we 
use the following shape tensor 



N 



Mjj = mkXj(k)xj(k), i,j = 1,2,3, 



(21) 



k=l 



where is the mass of the fc'th member particle and 
Xi is its position with respect to the halo center of mass. 
The shape tensor M is diagonalized for each halo. The 
ratios of principle axes are the square roots of ratios of 
the eigenvalues, and the orientations of the axes are the 
corresponding eigenvectors. However, the volume of the 
halo cannot be obtained from equation (|21|l . Conserv- 
ing the ellipsoidal shape of the halo, we scale up and 
down the shape to include almost all member particles. 
During the scaling, we set the virial volume to include 
95% of member particles since the halo boundary is not 
smooth. We check whether or not the virial region of a 
halo overlaps with those of other halos. If overlapping 
with each other, those halos are merged into a single 
halo. This merging process is similar to the linking of 
the FoF method. A halo found by this FoF-like method 
is a spatially isolated system whose virial region is not 
shared with other halos. We call this halo an isolated 
halo. 

In Figure |H1 we show the halo mass functions of the 
ACDM model simulation at z = measured from the 
distinguished halo (open circles) and the isolated halo 
catalogs (crosses) . Also shown are other analytic and fit- 
ting functions, such as the Press & Schechter (1974, solid 
curve), Sheth & Tormen (1999, dotted curve), Lee & 
Shandarin (1998, short-dashed curve), and Jenkins et al. 
(2001, long-dashed curve) mass functions. In the mass 
range of M h < 2 x 10 14 /i _1 M Q , the FoF-like method un- 
derestimates the number of halos, and overestimates the 
number of halos with mass Mh > 2 x 1O 14 /i _1 M0 when 
compared with the PSB catalog. Small distinguished ha- 
los embedded in dense regions are merged into large ha- 
los and, as a result, the mass function of isolated ha- 
los approximately follows the predicted mass function of 
Jenkins et al. (2001) who used the FoF method to find 
isolated halos in the Hubble Volume Simulations (Col- 
berg et al. 2000). We will study the this isolated and 
distinguished halos in the forthcoming paper in more de- 
tails. 

8. CONCLUSIONS 

We provide a new halo finding algorithm that effi- 
ciently and accurately identifies halos located in dense 
environments. This method employs a tidal radius con- 
straint as well as the total energy criterion in determining 
the halo membership . Based on the density field on a 
fine mesh with cell size of l\ = 2e, we detect halos with 
core size down to the force resolution e. We have com- 
pared the PSB method with other popular halo finding 
methods in resolving small halos embedded in a larger 
halo. The FoF and SO methods fail to identify the small 
halos. The tidal radius of a halo is a physically meaning- 
ful boundary descriptor, particularly in crowded regions. 
Density-related descriptors are not enough to identify ha- 
los when they are overlapping each other. When two 
halos are close to each other, the halos found by the 
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Fig. 8. — Halo mass functions obtained from the ACDM model simulation at z = 0. The halo mass functions obtained by the PSB 
(open circles) and the FoF-like method (crosses) are shown. We also plot the analytic and fitting functions proposed by Press & Schechter 
(solid curve), and Sheth & Tormen (dotted curve), Lee & Shandarin (short-dashed curve), and Jenkins et al. (long-dashed curve). 



HOP and DENMAX methods may have wrong physical 
characteristics such as the total mass, shape, spin, pe- 
culiar velocity and etc. The adaptive smoothing of the 
HOP method is not consistent with N-body simulations 
since the gravitational smoothing length e is kept con- 
stant throughout whole simulation. It is more reasonable 
to have the smoothing length h of the interpolating ker- 
nel W4 equal to e. The PSB algorithm is faster than the 
HOP because of the static smoothing length in obtaining 
the density field. 

The PSB method has several parameters. The most 
important ones are 5loc and S p . We use the constraint 
of <5loc <^ $vir < S p to find virialized halos with peak 
density greater than S p and to accurately calculate the 
tidal effect in the local overdense regions (S > ^loc)- 
With the parameter choices of (<5 p ,<5loc) = (312.5,3- 
5(72(0)), (312.5,(72(0)) for the concordance ACDM and 
SCDM models, we can identify halos with more than 30 
and 58 member particles in each model, respectively. We 
use the constraint of N core = 10 to avoid picking up ac- 
cidentally clustered unphysical peaks. There is another 
parameter Ni eve i , the number of density levels which split 
particles of a local group into at least Ni eve i particle 
sets. It is desirable for each particle set to be as small as 
possible to reduce underestimation of the halo potential 



(see Eq. We have found that there is no signifi- 

cant difference in the number of halo member particles 
if Ni eve i > 4. A parameter dependence study shows that 
the halo mass function obtained by the PSB method de- 
pends sensitively only on S p at the low mass end. Other 
parameters such as <5loCj Ni eve i, and N core7 are rather 
methodological parameters and basically do not affect 
the halo finding results. 

The main drawback of the PSB algorithm is its mem- 
ory requirement for the fine mesh. We need two fine 
meshes; one is for the density map and the other is for 
the linked list of particles. In 32-bit machines, the max- 
imum allowable number of fine cells is about 750 3 which 
need 3.4 GB memory budget. During the halo finding 
in the 1-Giga particle simulations, we have found that a 
few local-cell groups require more memory than the 32- 
bit limitation. In principle, there is no memory limita- 
tion problem on 64-bit machines as far as there is enough 
physical memory allocated to each CPU. The efficiency 
of the PSB algorithm decreases when the number of ha- 
los increases in a local particle group. The increasing 
number of halos induces the computational overload to 
measure the tidal radius with respect to all other mas- 
sive halos. Therefore, the measurement of a tidal radius 
scales as 0(N%). We find that about 30% of the CPU 
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time in the halo finding is spent in this job. Algorithms 
faster than the current version of the PSB algorithm are 
desirable for the next largest simulations. 

We apply our method to the 1-Giga particle simula- 
tions of the ACDM and SCDM models. About four hours 
are taken for the halo finding in the simulations of the 
ACDM and SCDM models at z = on 32 IBM POWER- 
4 chips. We expect the PSB method easily scalable to 
2048 3 or even more massive N-body simulations due to 
its algorithm of dividing whole simulation particles into 
multiple local particle groups. 
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